﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using DotSpatial.Data;
using DotSpatial.Projections;
using DotSpatial.Topology;
using System.Threading.Tasks;

namespace GeostatisticalTool.Classes
{
    public class SpatialContinuity : ISpatialContinuityData
    {
        private DataSourceInterpolators dataMethod;
        private int npoints = 0;
        private List<Kpoint> kpoints;
       // private double maux = 0;
       // private double mauxV = 0;
        private Lag lag=null;
        private Stat[,] spatConti;
        private Stat[,] spatDist;
        private Covariance[,] spatContiCov;
        private Stat semivariance;
        private Stat distance;
        private Stat covariance;
        private SortedList<int, double> distancesA = null;
        private SortedList<int, double> semivarianceA = null;
        private SortedList<int, double> covarianceA = null;
        private bool calculated = false;

       // public event EventHandler ChangedSemivariances;
        public event EventHandler ChangedLagParameter;
       
        
        protected virtual void OnChangedLagParameter(EventArgs e)
        {
            if (ChangedLagParameter != null)
                ChangedLagParameter(this, e);

            calculated = false;
        }

        


        #region properties

        public Stat[,] DistancesMeanOfGrid
        {
            get { return spatDist; }
        }

        public Covariance[,] CovarianceOfGrid
        {
            get { return spatContiCov; }
        }

        public Stat[,] SemivarianceOfGrid
        {
            get { return spatConti; }
        }


        public Stat SemivarianceStatistic
        {
            get { return this.semivariance; }

        }

         public Stat DistanceStatistic
        {
            get { return this.distance; }
        }

         public Stat GetDistanceStatistic()
         {
             return this.distance;
         }

         public Stat CovatiranceStatistic
        {
            get { return this.covariance; }
        }

        public Lag LagValue
        {
            get
            {
                return lag;
            }

            set
            {
                if (this.LagSize != value.lagSize ||
                    this.NumberOfLags != value.numLags ||
                    this.lag.tolerance != value.tolerance)
                {
                    this.lag = value;
                    Grid();
                    CalculateDistancesLag();
                    calculated = false;
                    OnChangedLagParameter(EventArgs.Empty);
                }
            }


        }

        public int NumberOfLags
        {
            set
            {
                this.NumberOfLags = value;
                Grid();
                CalculateDistancesLag();
                OnChangedLagParameter(EventArgs.Empty);
            }
            get
            {
                return this.lag.numLags;
            }
        }

        public double LagSize
        {
            set
            {
                this.LagSize = value;
                Grid();
                CalculateDistancesLag();
                OnChangedLagParameter(EventArgs.Empty);
            }

            get
            {
                return this.lag.lagSize;
            }
        }

        public double[] GetSemivariances
        {
            get
            {
                if (calculated)
                    return semivarianceA.Values.ToArray();
                else
                {
                    CalculateListValues();
                    return semivarianceA.Values.ToArray(); 
                }
            }
        }

        public double[] GetDistances
        {
            get
            {
                if (calculated)
                    return distancesA.Values.ToArray();
                else
                {
                    CalculateListValues();
                    return distancesA.Values.ToArray(); ;
                }
            }
        }

        public double[] GetCovariances
        {
            get
            {
                if (calculated)
                    return covarianceA.Values.ToArray();
                else
                {
                    CalculateListValues();
                    return covarianceA.Values.ToArray();
                }
            }
        }
        
        #endregion


        public SpatialContinuity(DataSourceInterpolators data,Lag lag)
        {

            this.semivariance = new Stat(true);
            distance = new Stat(true);
            covariance = new Stat(true);

            this.dataMethod = data;
            kpoints = data.Kpoints;
            npoints = kpoints.Count;
            this.lag = lag;
            Grid();
            CalculateDistances();
            CalculateDistancesLag();
            calculated = false;
        }


        public SpatialContinuity(DataSourceInterpolators data)
        {

            this.semivariance = new Stat(true);
            distance = new Stat(true);
            covariance = new Stat(true);

            this.dataMethod = data;
            kpoints = data.Kpoints;
            npoints = kpoints.Count;
            CalculateDistances();
            this.lag = new Lag(12, distance.Max*0.7/12, 0.5);
            Grid();
            CalculateDistancesLag();
            calculated = false;
        }


        private void Grid() {

            if (lag == null)
                return;


           spatConti= new Stat[lag.numLags*2,lag.numLags*2] ;
           spatContiCov = new Covariance[lag.numLags * 2, lag.numLags * 2];
           spatDist = new Stat[lag.numLags * 2, lag.numLags * 2];
            for(int row=0; row<lag.numLags*2;row++ )
            {
              for(int col=0; col<lag.numLags*2;col++ )
              {
                 spatConti[row,col]=new Stat(true);
                 spatContiCov[row, col] = new Covariance(true);
                 spatDist[row, col] = new Stat(true);
              }
            }
        }

        public IFeatureSet InfluenceZone(string majorR, string minorR, string az)
        {  
            double value=1;
            try{
               value=Convert.ToDouble(majorR);
               if (!minorR.Equals(""))
               { 
                   return  InfluenceZone(value,Convert.ToDouble(minorR), Utils.ToRad(Convert.ToDouble(az)));
               }
            } catch
            {
            
            }


        return InfluenceZone(lag.numLags * lag.lagSize,lag.numLags * lag.lagSize,0);
        
        }

        public IFeatureSet InfluenceZone(double majorR, double minorR, double az)
        {
            if (lag == null)
                return null;

            FeatureSet border = new FeatureSet(DotSpatial.Topology.FeatureType.Polygon);
            border.Projection = KnownCoordinateSystems.Projected.UtmWgs1984.WGS1984UTMZone10N;
            double size = lag.numLags * 2 * lag.lagSize;

            Coordinate center = new Coordinate(500000 + (size / 2), size / 2);

            //Coordinate[] coorPol = new Coordinate[73];
            //int j=0;
            //for (int i = 0; i < 360; i += 5)
            //{ 
            
            //coorPol[j]=AzimuthDist(center, i*Math.PI/180.0,lag.numLags  * lag.lagSize);
            //    j++;
            //}

            //coorPol[72] = AzimuthDist(center, 0, lag.numLags * lag.lagSize);
            //LinearRing pol = new LinearRing(coorPol);
            //Polygon poly = new Polygon(pol);
            border.Features.Add(Utils.DrawEllipse(new Kpoint(500000 + (size / 2), size / 2,0,0),majorR,minorR, az));

            
            return border;
        }


        private void CalculateListValues()
        {
            if (lag == null)
                return;

            distancesA = new SortedList<int, double>();
            semivarianceA = new SortedList<int, double>();
            covarianceA = new SortedList<int, double>();
            double size = lag.numLags * 2 * lag.lagSize;

            int nr = (lag.numLags * 2) ;
            int nc = (lag.numLags * 2) ;
            int i = 0;
            for (int row = 0; row < nr; row++)
            {
                for (int col = 0; col < nc; col++)
                {

                   if (col >= (lag.numLags * 2) || row >= (lag.numLags * 2))
                    { }
                    else
                    {
                        if (spatConti[row, col].num > 0 && GetDistanceGrid(row, col) <= (lag.numLags  * lag.lagSize))
                        {

                        // distancesA.Add(i, spatDist[row, col].Mean);
                        distancesA.Add(i, GetDistanceGrid(row, col));
                        semivarianceA.Add(i, spatConti[row, col].Mean);
                        covarianceA.Add(i, spatContiCov[row, col].Result());
                        i++;
                    }
                }
            }
            calculated = true;
        }

   }

        public IRaster Surface( string outGridPath)
        {
            if (lag == null)
                return null;

            distancesA = new SortedList<int, double>();
           
            semivarianceA = new SortedList<int, double>();
             covarianceA = new SortedList<int, double>();
            double size = lag.numLags * 2 * lag.lagSize;
            Extent ext = new Extent(new double[4] {500000,0,500000+size,size });
            IRaster raster;
            AbstractInterpolator.CreateGridFromExtents(ext, lag.lagSize, KnownCoordinateSystems.Projected.UtmWgs1984.WGS1984UTMZone10N, -1, outGridPath, out raster);

            int nr = raster.NumRows;
            int nc = raster.NumColumns;
            int i = 0;
            for (int row = 0; row < nr; row++)
            {
                for (int col = 0; col < nc; col++)
                {
                    if (col >= (lag.numLags * 2) || row >= (lag.numLags * 2))
                    { }
                    else
                    {
                        if (spatConti[row, col].num > 0 && GetDistanceGrid(row, col) <= (lag.numLags  * lag.lagSize))
                        {

                            raster.Value[row, col] = spatConti[row, col].Mean;
                           // distancesA.Add(i, spatDist[row, col].Mean);
                            distancesA.Add(i,GetDistanceGrid(row,col));
                            semivarianceA.Add(i, spatConti[row, col].Mean);
                            covarianceA.Add(i, spatContiCov[row, col].Result());
                            i++;
                        }
                    }
                }
            }
            raster.Save();
            raster.Close();
            calculated = true;
            return raster;
        }

        public double GetDistanceGrid(int row, int col)
        {
            if (lag == null)
                return 0;

            double initialcol = this.LagSize / 2;
            double initialrow = this.LagSize / 2;
            double len = this.NumberOfLags * this.LagSize;
            double x,y;

                x = (len) - initialcol - ((col ) * this.LagSize);
                y = (len) - initialcol - ((row) * this.LagSize);

                return Math.Sqrt((x * x) + (y * y));
        }


        public Coordinate AzimuthDist(Coordinate ptoI, double azi, double dist)
        {
                return new Coordinate(ptoI.X + Math.Sin(azi) * dist, ptoI.Y + Math.Cos(azi) * dist);
        }

       

        private double GetAzimuth(int i, int j)
        {
            //double dx = (this.cpoints[j, 0] - this.cpoints[i, 0]);
            //double dy = (this.cpoints[j, 1] - this.cpoints[i, 1]);
            double dx = this.kpoints[j].X - this.kpoints[i].X;
            double dy = this.kpoints[j].Y - this.kpoints[i].Y;
            if (dx == 0 && dy > 0) return 0;
            if (dx == 0 && dy < 0) return (Math.PI);
            if (dx > 0 && dy == 0) return (Math.PI / 2);
            if (dx < 0 && dy == 0) return (3 * Math.PI / 2);
            if (dx > 0 && dy > 0) return (Math.Atan(dx / dy));
            if (dx > 0 && dy < 0) return (Math.PI + Math.Atan(dx / dy));
            if (dx < 0 && dy < 0) return ((Math.PI) + Math.Atan(dx / dy));
            if (dx < 0 && dy > 0) return ((2 * Math.PI) + Math.Atan(dx / dy));
            return 0;
        }

        /// <summary>
        /// Calculate the distastance between two points
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>Return the distance</returns>
        private double Mdistances(
            int i, int j)
        {
            return Math.Sqrt(((this.kpoints[i].X - this.kpoints[j].X) * (this.kpoints[i].X - this.kpoints[j].X)) +
                             ((this.kpoints[i].Y - this.kpoints[j].Y) * (this.kpoints[i].Y - this.kpoints[j].Y)));
        }


        /// <summary>
        /// Calculate the distastance between two points
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>Return the distance</returns>
        private double[] MdistancesXY(
            int i, int j)
        {
            return new double[2] {(this.kpoints[j].X - this.kpoints[i].X) ,  (this.kpoints[j].Y - this.kpoints[i].Y)};
        }

        /// <summary>
        /// Calculate Semivariances
        /// </summary>
        /// <param name="i">Position ith in the point array</param>
        /// <param name="j">Position jth in the point array</param>
        /// <returns>the semivariance</returns>
        private double MsemiVarEmp(
             int i, int j)
        {
            return 0.5 * ((this.kpoints[i].W - this.kpoints[j].W) * (this.kpoints[i].W - this.kpoints[j].W));
        }

         /// <summary>
        /// This method calculate the distances between all points
        /// </summary>
        private void CalculateDistances()
        {
            double[] mdisv;
            double mdis;
            double msem = 0d;
           Parallel.For(0,this.npoints, i =>
           {
               Parallel.For(0, this.npoints, j =>
             {
                 mdisv = this.MdistancesXY(i, j);
                 msem = this.MsemiVarEmp(i, j);
                 mdis = this.Mdistances(i, j);

                 if (mdisv[0] != 0 && mdisv[1] != 0)
                 {
                     Stat semi = new Stat(msem);
                     Covariance cov = new Covariance(this.kpoints[j].W, this.kpoints[i].W);
                     Stat dist = new Stat(mdis);

                     distance += dist;
                     this.semivariance += semi;

                 }
             });
			});
			

            //for (int i = 0; i <this.npoints ; i++)
            //{
            //    for (int j = i; j < this.npoints; j++)
            //    {
            //        mdisv = this.MdistancesXY(i, j);
            //        msem = this.MsemiVarEmp(i, j);
            //        mdis = this.Mdistances(i, j);

            //        if (mdisv[0] != 0 && mdisv[1] != 0)
            //        {
            //            Stat semi = new Stat(msem);
            //            Covariance cov = new Covariance(this.kpoints[j].Z, this.kpoints[i].Z);
            //            Stat dist = new Stat(mdis);

            //            distance += dist;
            //            this.semivariance += semi;

            //        }
            //    }
            //}
        
        }




        /// <summary>
        /// This method calculate the distances between all points
        /// </summary>
        private void CalculateDistancesLag()
        {
            double[] mdisv ;
            double mdis;
            double msem = 0d;

            
            Parallel.For(0,this.npoints, i =>
           {

               Parallel.For(0, this.npoints, j =>
              {
                  mdisv = this.MdistancesXY(i, j);
                  msem = this.MsemiVarEmp(i, j);
                  mdis = this.Mdistances(i, j);

                  if (mdisv[0] != 0 && mdisv[1] != 0)
                  {


                      int[] cells = lag.Cells(i, j, mdisv[0], mdisv[1], GetAzimuth(i, j));

                      int lim = lag.numLags;
                      Stat semi = new Stat(msem);
                      Covariance cov = new Covariance(this.kpoints[j].W, this.kpoints[i].W);
                      Stat dist = new Stat(mdis);

                      //distance += dist;
                      //semivariance += semi;
                      //  covariance += new cov;
                      int v = (this.NumberOfLags * 2) - 1;
                      if (CheckBorder(cells[2]) && CheckBorder(cells[0]))
                      {
                          spatConti[cells[2], cells[0]] += semi;
                          spatContiCov[cells[2], cells[0]] += cov;
                          spatDist[cells[2], cells[0]] += dist;

                          spatConti[v - cells[2], v - cells[0]] += semi;
                          spatContiCov[v - cells[2], v - cells[0]] += cov;
                          spatDist[v - cells[2], v - cells[0]] += dist;

                      }

                      if (CheckBorder(cells[3]) && CheckBorder(cells[0]))
                      {
                          spatConti[cells[3], cells[0]] += semi;
                          spatContiCov[cells[3], cells[0]] += cov;
                          spatDist[cells[3], cells[0]] += dist;

                          spatConti[v - cells[3], v - cells[0]] += semi;
                          spatContiCov[v - cells[3], v - cells[0]] += cov;
                          spatDist[v - cells[3], v - cells[0]] += dist;

                      }
                      if (CheckBorder(cells[2]) && CheckBorder(cells[1]))
                      {
                          spatConti[cells[2], cells[1]] += semi;
                          spatContiCov[cells[2], cells[1]] += cov;
                          spatDist[cells[2], cells[1]] += dist;

                          spatConti[v - cells[2], v - cells[1]] += semi;
                          spatContiCov[v - cells[2], v - cells[1]] += cov;
                          spatDist[v - cells[2], v - cells[1]] += dist;
                      }
                      if (CheckBorder(cells[3]) && CheckBorder(cells[1]))
                      {
                          spatConti[cells[3], cells[1]] += semi;
                          spatContiCov[cells[3], cells[1]] += cov;
                          spatDist[cells[3], cells[1]] += dist;

                          spatConti[v - cells[3], v - cells[1]] += semi;
                          spatContiCov[v - cells[3], v - cells[1]] += cov;
                          spatDist[v - cells[3], v - cells[1]] += dist;
                      }
                  }

              });
            });
        }

        private bool CheckBorder(int v)
        {
            if (v != int.MinValue)
            {
                if (v >= 0 && v <=((lag.numLags*2)-1))
                    return true;
                else
                    return false;
            }
            else
            {
                return false;
            }
        }


        public double[] GetDistancesValues()
        {
            return this.GetDistances;
        }

        public double[] GetSemivariancesValues()
        {
            return this.GetSemivariances;
        }


        public Stat GetSemivarianceStatistic()
        {
            return this.SemivarianceStatistic;
        }

        public Lag GetLag()
        {
            return this.LagValue;
        }
    }


}
